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Abstract 

We present a scheme to solve the nonlinear multigroup radiation dif- 
fusion (MGD) equations. The method is incorporated into a massively 
parallel, multidimensional, Eulerian radiation-hydrodynamic code with 
adaptive mesh refinement (AMR). The patch- based AMR algorithm re- 
fines in both space and time creating a hierarchy of levels, coarsest to 
finest. The physics modules are time-advanced using operator splitting. 
On each level, separate "level-solve" packages advance the modules. Our 
multigroup level-solve adapts an implicit procedure which leads to a two- 
step iterative scheme that alternates between elliptic solves for each group 
with intra-cell group coupling. For robustness, we introduce pseudo tran- 
sient continuation ($tc). We analyze the magnitude of the ^Ptc parameter 
to ensure positivity of the resulting linear system, diagonal dominance and 
convergence of the two-step scheme. For AMR, a level defines a subdo- 
main for refinement. For diffusive processes such as MGD, the refined level 
uses Dirichet boundary data at the coarse-fine interface and the data is de- 
rived from the coarse level solution. After advancing on the fine level, an 
additional procedure, the sync-solve (SS), is required in order to enforce 
conservation. The MGD SS reduces to an elliptic solve on a combined 
grid for a system of G equations, where G is the number of groups. We 
adapt the "partial temperature" scheme for the SS; hence, we reuse the 
infrastructure developed for scalar equations. Results are presented. We 
consider a multigroup test problem with a known analytic solution. We 
demonstrate utility of ^tc by running with increasingly larger timesteps. 
Lastly, we simulate the sudden release of energy Y inside an Al sphere 
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(r = 15 cm) suspended in air at STP. For Y = 11 kT, we find that gray 
radiation diffusion and MGD produce similar results. However, if Y = 1 
MT, the two packages yield different results. Our large Y simulation con- 
tradicts a long-standing theory and demonstrates the inadequacy of gray 
diffusion. 



1 Introduction 

This paper describes a numerical method to solve the radiation multigroup 
diffusion (MGD) equations. Two themes are presented. One is the scheme 
itself. We add Pseudo Transient Continuation (Wtc) to the familiar "fully im- 
plicit" method of Axelrod et al [2]. The second theme is code-specific. Our 
MGD solver is embedded in a multidimensional, massively parallel, Eulerian 
radiation-hydrodynamic code, which has patch-based, time-and-space Adaptive 
Mesh Refinement (AMR) capability. Our code's AMR framework stems from 
the Berger and Oliger idea [3] developed for hyperbolic, compressible hydrody- 
namic schemes. The idea was expanded by Almgren et al [1] and applied to 
the type of elliptic solvers required for the incompressible equations of Navier- 
Stokes. Howell and Grecnough [7] applied the Almgren et al framework to the 
scalar, parabolic "gray" radiation diffusion equation, thereby creating the start 
of our radiation-hydrodynamic code. 

The AMR framework works as follows. A domain, referred to as the "coarse" 
or LO level, is discretized using a uniform, coarse spatial mesh size h c } After 
advancing with a timestep At c , the result is scanned for possible improvement. 
One may refine subregions containing a chosen material, at material interface(s), 
or at shocks, etc. Whatever refinement criteria are used, after the subdomains 
are identified, specific routines define a collection of "patches," which cover the 
subdomains. In two dimensions, the patches arc unions of rectangles; in 3D, 
they are unions of hexahedra. The patches need not be connected, but they 
must be contained within the coarse level. The patches denote the "fine" or 
LI level and are discretized with a uniform, spatial mesh size hf. A typical 
refinement ratio h c /hf equals two, but higher multiples of two are also allowed. 

Because the original framework was designed for temporally explicit hyper- 
bolic schemes, At c is restricted by a CFL condition. This implies a similar 
restriction for the LI level timestep Atf. For the case, h c /hf — 2, level LI 
time-advances twice using Atf = At c /2. Boundary conditions for level LI are 
supplied as follows. Wherever level LI extends to the physical boundary, the 
level uses the conditions prescribed by the problem. Portions of level Li's bound- 
ary which lie inside the physical domain have conditions prescribed by time and 
space interpolated data obtained from the LO solution. For diffusion equations, 
these conditions are of Dirichlct type. The numerical solution consists of both 
coarse and fine grid results. Unfortunately, as it stands, the composite solution 
does not guarantee conservative fluxes across the level boundaries. To maintain 



1 In multiple dimensions, coordinates have their own mesh spacing. 
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conservation, a separate procedure, dubbed a sync-solve (SS) is required. The 
SS reduces to an elliptic unstructured grid solve on the composite grid of LO and 
LI levels. The AMR procedure may be recursive. That is, a level LI grid may 
generate its own subdomain for refinement, i.e., a level L2. In that case, one SS 
couples results from levels LI and L2. Once the levels advance to the LO level 
time, a SS coupling all three levels ensues. For the multigroup equations, the 
SS requires an unstructured grid solve for a coupled system of reaction-diffusion 
equations. Our scheme for a multigroup SS is an important theme of this paper. 

The MGD equations stem from a discretization of the multifrcqucncy ra- 
diation diffusion equations. The latter is an approximation to the equations 
of radiation transfer, obtained by assuming the matter to be optically thick, 
which suppresses the directional dependence of the radiation intensity. Details 
of the derivation may be found in various sources: Mihalas and Mihalas [12], 
Zel'dovich and Raizer [26], Pomraning [17]. 

The gray radiation diffusion equation is a simplification of the MGD equa- 
tions. It is essentially a one-group equation and is derived by integrating over all 
frequencies. Surprisingly, it gives very good results in many cases. However, it 
clearly cannot display frequency-dependent effects. When those are important, 
it gives incorrect results. Unfortunately, unless one solves a problem with both 
gray and MGD, one never knows when the former is adequate. 

We now summarize the paper. Our MGD scheme consists of two parts. 
Sections 2 and 3 develop the level-solve algorithm, which is applied on each 
level. Section 2 develops the equations, the discretization, and our \Ptc scheme. 
Section 3 proves three lemmas which determine the initial magnitude of the Vl^tc 
parameter a. Our philosophy for a is as follows. The result of the level solve is 
the time-advanced radiation group energy density, which physics dictates to be 
nonnegative. Zeroing anomalously negative values is not an option since they 
are the correct conservative solution to the linear system that stems from the 
discretization of the system. Thus, the unphysical result nonetheless conserves 
energy. The difficulty is avoided if in the original formulation of the linear system 
Ax = b, A is an M-matrix and the right-hand-side (RS) is nonnegative. Since 
we solve Ax = b using an iterative scheme, the magnitude of a is determined 
to ensure b > 0, a diagonally dominant A, and that the iterations converge. 
To a large extent, we are guided by Pert [16], who discusses how and why the 
solution to a discretization of an equation may be unacceptable from a physical 
standpoint. For a first reading, section 3 may be skipped; the analysis of the 
required magnitude of a is not needed for the subsequent sections. 

We note that \ftc is widely used to solve nonlinear systems of equations. 
It is closely related to the Inexact Newton Backtracking Method by Shahid et 
al [18]. When applying "ftc to a Newton solver, the basic idea is to limit the 
change to the iterates when one is far from the root but not restrict the change 
as one approaches the root. With Wtc, limiting is done by the magnitude of the 
pseudo-timcstcp. Kcllcy and Keyes [8] put \I>tc on a solid analytic framework by 
examining the three regimes of "ftc: small, medium, and large pseudo-timesteps. 
In the last regime, ^>tc recovers Newton's second order of convergence. 

Our \Ptc implementation differs from the norm. Standard applications typ- 
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ically detect when a problem is "hard" and then reduce the timcstcp or some 
other parameter by an arbitrary amount. However, this method will not work 
for us because our solver is embedded in a time-dependent multiphysics code 
with separate modules for compressible gasdynamics, heat conduction, radia- 
tion transport. Our MGD solver is called numerous times during the course of a 
simulation. (If running, with AMR, multiple times per physical time advance.) 
Although the physical At is controlled by various means, and depending on the 
problem can vary many orders of magnitude, we require a MGD solver that 
works under all conditions. Our vl/tc approach is similar to the one of Shestakov 
et al [20]. We set the initial magnitude of the "ftc parameter to ensure that 
for the first step, our iteration scheme converges and that the result is physical. 
We note that our usage of \ttc is nearly equivalent to having the MGD module 
time-advance not in a single (physical) step At, but in smaller time increments 
until the desired time t° + At is reached. Some colleagues refer to the process 
as "sub-cycling" the radiation module. It is easy to show that the lemmas of 
Sec. 3 still apply for sub-cycling. 

Section 4 describes the second part of our solver, viz., the sync-solve. Sec- 
tion 5 contains results. Three problems are presented. The first, in Sec. 5.1, 
displays the accuracy of the method and its convergence properties: first order 
in time and second order in space. Section 5.2 demonstrates the utility afforded 
by SSftc. For hard problems, it accelerates convergence; for very hard problems, 
iPtc is indispensable. Section 5.3 models the explosive expansion of a hot metal 
sphere suspended in cold air. The simulation couples all of the code's physics 
modules. The problem is an ideal candidate for AMR since effects propagate 
a large distance away from the source, yet in early times, resolution is needed 
only near the sphere. The problem also demonstrates the necessity of multi- 
group diffusion. We find that if the sphere's energy is very high, gray diffusion 
gives the wrong answer. For a 1 MT energy source, our MGD simulation contra- 
dicts results of Brode [5], who used gray diffusion. Section 6 contains concluding 
remarks. 

There are three appendices. Appendix A gives a table of exact values for 
the test problem described in section 5.1. Appendix B discusses situations 
that may complicate attaining a diagonally dominant matrix when discretizing 
the multigroup system. Appendix C presents a spatial convergence analysis 
of the multigroup system when running in "production" mode, that is, with a 
dominant flux limiter and with AMR. 



Ignoring velocity terms and Compton scattering, the multifrcquency radiation 
equations (CGS units) (Mihalas and Mihalas [12]) are: 
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In (l)-(2), u v and e represent the spectral radiation energy density and matter 
specific energy respectively. The former is a function of position x, time t and 
frequency v, while e is a function of the mass density p and material temperature 
T, quantities which themselves depend on x and t. Evolution of p is governed 
by hydrodynamics. Hence, in our context, p is a known function. Introducing 
the specific heat c v — de/dT turns (2) into an evolutionary equation for T; 
hence, the left-hand-side (LS) becomes pc v dtT. The subscript v designates that 
the term varies with frequency. In (l)-(2), c denotes the speed of light, n v the 
absorption opacity, and B v the Planck function, 

B v = (8?r h/c 3 ) v 3 / [exp(hv/kT) - 1] (erg sec cm" 3 ) , 

where h and k are the Planck and Boltzmann constants, respectively. The 
diffusion coefficient D v depends on the total inverse mean free path \v = P^ v + 
pK. v ,s! where n v and k„ jS are the absorption and scattering opacities, respectively. 
(The opacities are also functions of material composition, p and T.) In (1), 
the term — D v Wu denotes the spectral radiation energy flux. To limit energy 
streaming faster than c, a flux limiter is introduced, e.g., 

A, = c/[3x„ + |V(u„)|/u„]. (3) 

The multigroup equations are derived as follows. The frequency domain is 
discretized into G groups with boundaries {f ff }^ =0 satisfying 

< Uq < V\ < . . . < v g < oo . 

The equations are integrated over groups. We define 

u g (x, t) = u v = dvu v . 

Jg Jvg-i 

Time derivatives are replaced by differences and the system is multiplied by the 
timestcp At. Integration of the transport and absorption terms requires defin- 
ing group-averaged opacities. Linearizing the Planck function about a known 
temperature T* , the absorption term is expressed as 

k„ (B v - u„) = K g [B g + B' g (T — T*) — u g ], 

where n g is the group-averaged absorption opacity, B g — J g B u \t=t* , and 
B' g = f g {dB v /dT)\ T = T ,. In a semi-implicit scheme, T* = T°, where T° is 
the temperature at the start of the time cycle. For fully implicit differencing, 
we must iterate until T* converges to T. For the transport term, we define 

At / V • D v Vm = V • D g \7u g , 

where D g depends on a group-averaged inverse mean free path \g- Note that 
At has been absorbed into D g . 
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The above definitions yield the multigroup equations, 

= U g — Ug + V ■ DgVUg + K g { u, T) , g=l,...,G (4) 
G 

= pc v (T -T)-^K e (u,T) , (5) 
1=1 

where u g and T° denote values at the start of the time-advance, 

K g (u,T) = a g [B g + B' g (T-T*)-u g ], 
a g = AtcpKg . 

Equations (4)-(5) comprise a nonlinear system with the strongest nonlin- 
earity due to the emission term B. To a lesser extent, opacities also have a 
temperature dependence and for nonideal gases, so does c v . However, for ease of 
solution, we may choose to view (4)-(5) as a linear system in which case all coeffi- 
cients are evaluated at the old temperature T° . For simulations in which matter 
and radiation are tightly coupled, i.e., where we expect to have u v = B v , the 
solution to the semi-implicit difference equations is u g = B g + B' g (T — T ) , with 
B g and B' g evaluated at T = T . For high frequencies, lim„_> 00 (.B„ / 'B' v ) ~ 1/f; 
hence, B g <C B' g for large g. Unfortunately, if the temperature is decreasing, 
i.e., if (T — T°) < 0, the linearized emission term is negative for large g, leading 
to the unphysical result: u g < 0. On the other hand, if we are able to iterate 
on T* so that it converges to T, then in tightly coupled simulations, we obtain 
the desired solution u g = B g with B g evaluated at the advanced temperature. 

In our code we provide both options, i.e., solving a linear system, or con- 
verging on the implicit source. 2 In either case, solving (4)-(5) on a large domain 
with many groups presents a formidable task. To facilitate the task, we intro- 
duce pseudo transient continuation (iPtc) and replace the zeros on the LS of 
(4)-(5) with the ^tc derivatives, 

T{ug-u* g ) and P c v t{T-T*) , 

where r > 0, the inverse of the pseudo-timestep, is the ^tc parameter whose 
magnitude is at our disposal. 

The variables u* and T* represent advances in pseudo time; they always 
appear on the LS of (4)-(5). As mentioned above, we provide the option of 
running in either semi-implicit (SI) or fully-implicit (FI) mode. With SI, since 
B g is linearized about T = T°, in the definition of the coupling term K g , we 
substitute T° for T* . However for FI, K g is defined as above; B g is linearized 
about the pseudo time temperature T* . The two modes lead to subtle differences 
in the scheme, as shown below. 

For the FI scheme, if the matter equation is solved for the temperature 
change, we obtain 

G 

( T - T* ) - pc„ (T° - t*) - (B e - u t ) , (6) 



2 At the time of this writing, opacities and c v were time-lagged. 
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where 

G 

= pc v a + ai B' g and a = 1 + t . (7) 

e=i 

The domain of relevance t > corresponds to cr > 1 . 
For the SI scheme, the temperature change is, 

G 

8- 1 (T-T°)=pc v (<j-l) (T* - T°) - £ a, (5/ - «*) . (8) 

The term S is defined as above, but B g and B' are evaluated at T = T°. 

For the FI scheme, if (6) is substituted into the equation for u g , we obtain 

G 

-V • D g Vu g + {a + a g )u g - fg^aeUj = 

u g + (a-l)u*g + a g Bg + fg(^pc v (T -T*)-]T>.B^ , (9) 

where f g = 8 a g B' g . Equation (7) implies f g < 1, for all 5. For the SI scheme, 
the RS of (9) changes: pc v (T° — T*) is replaced with pc„ (cr - 1) (T* - T°). 
Equation (9) corresponds to a linear system 

Au = w 

of order (TV x G), where N is the number of mesh cells and G the number 
of groups. The first term on the LS of (9) consists of second order, central 
differences over space. We write this term as 

-V • Dg VUg = +V d ,g Ug - V^g Ug . 

The first part represents multiplication of the vector u g by a diagonal matrix; 
the second term denotes multiplication by the off-diagonal part. The coefficients 
of T>d and V Q are nonnegative. 

On the LS of (9), the term —f g J2i=i a t u t IS referred to as the "re-emission 
source" [13], since it represents radiation energy absorbed by matter and re- 
emitted. If we define the column vectors / and a with components f g and a g , 
respectively, the re-emission term is expressed as the matrix- vector product 

-(fa T )u, (10) 

where a T = transpose (a), and u is the column vector of unknowns. Since the 
re-emission term does not couple cells, (10) corresponds to separate products: 
one per cell, with each product of order G. 

These observations allow expressing the matrix as 



A = A - Mi - M 2 , 



(11) 
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where A is diagonal, Mi contains the offdiagonal terms due to the (spatial) 
diffusion term, and M 2 contains the offdiagonal terms due to intcrfrcquency 
coupling. The corresponding elements are 

Ag = V^g + a + Clg - fgClg , 

(Ml U)g = V ,g Ug , 

G 

(M 2 u) g = fg^2a e u e . 

The decomposition (11) leads to the iterative scheme proposed by Axelrod et al 
[2], which improves a guess by successively solving 

(A- M 2 )u {t+1 'V = «; + MiuW (12) 
(A-Mi)u (i+1) = w + M 2 u {t+1/2) . (13) 

We solve (12)-(13) until t*W converges. Convergence is gauged by evaluating 
the 1-norms of w and the residual r = w — Au; the latter defined as, 

r = w- Au^+V = M 2 (u^+V - U (*+V2)) . 

The procedure is fast since multiplication by M 2 is local to each cell, which is 
very convenient if the spatial domain is decomposed on multiple processors. 

We now review the derivation of the system Au = w. First, we assume that 
^tc is not used, i.e., that a = 1 in (7)-(9). For the SI scheme, the terms Bi 
and B' e are evaluated at T = T°. For FI differencing, we require two types of 
iterations. Equations (12) and (13) comprise the inner iteration. It is initialized 
with z/ ) equal to u°. Once the inner iteration has converged to sufficient accu- 
racy, (6) yields the new temperature. The SI scheme essentially ends after the 
inner iteration converges (see below). For FI differencing, after T is computed, 
the outer iteration sets T* = T, recomputes Bg and B' e at T — T* and returns 
to the inner iteration. The outer iteration halts when T* converges. 

If ^Ptc is invoked, more care is required because when a > 1, the system 
Au = w is not a true discretization of the multigroup equations. Despite this 
complication, \ftc brings robustness to the scheme. The ^tc parameter r plays 
the role of an inverse timestep in pseudo-time. In principle, we could set t to a 
large value and solve a succession of linear systems. The solution of each system 
represents an advance in pseudo-time. We continue advancing until we reach 
the pseudo-time steady-state. This is easily seen by letting u* = u g on the RS of 
(9) and moving the term to the LS. However, making r large is not practical as 
it involves many pseudo-time advances. Furthermore, the intermediate pseudo- 
time results are of no interest. Consequently, we adopt the strategy of making 
r as small as possible. We discuss the strategy in section 3. 

^tc may be used with either SI or FI differencing. In the former, once (12) 
and (13) are converged, (8) yields the new temperature T. We then compute 
the 1-norm of the "nonlinear" residual of the linearized equation for the matter 
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energy, 

rm = v(^pc v (T- T°) - ^ a e [B t + B' e (T- T°) - u,]^ , (14) 

where U is the cell volume. The residual is compared with the 1-norm of the 
matter "energy" V pc v T, and in order to monitor stagnation, it is also compared 
with the energy change over the pseudo-timestep Vpc v (T — T*). With FI, the 
temperature T, obtained from (6), is used to compute the emission Bg. The 
residual r n i is defined as in (14), except without the B' t (T - T°) term. 

Unfortunately, unless the iterations converge to round-off accuracy, energy 
may not be conserved. Lack of conservation stems from values of user-set 
parameters that control stopping criteria for the iterations. For example, it 
may be efficient to halt once ||r n j||i < 1CP 6 , and the norm of the iterates 
||(AT)/T|| 00 < 10~ 2 since continuing brings little noticeable (visual) improve- 
ment to the solution. However, if one were to stop at that point, energy may 
not be conserved to desired accuracy. To restore conservation, we provide the 
option of an additional step. After the iterations stop, we assume that the last 
computed temperature T is "frozen" and use it to compute emission. In the 
SI scheme, emission into the gth group is defined as S g = B g + B'(T — T°), 
where B g and B' g are evaluated using T°. (To prevent unphysical behavior, S g 
is not allowed to be negative.) In the FI scheme, we evaluate B g using T and set 
S g = B g . Having a known emission allows us to compute the energy-conserving 
radiation field. The groups decouple. For g = 1, . . . , G, we solve 

-V • D g Vug + {1 + a g )u g = u° g + a g Sg . 

After computing u g , the matter energy density change is 

G 

AS = -^ai{B g - u g ) , 

i=i 

where, if using the SI scheme, B g is linearized about T = T°, or with FI, is 
evaluated at T. The quantity V AS represents the average energy change of 
the matter. In cells with more than one material, we adapt a suggestion of 
Zimmerman [28], which simulates intra-cell gray diffusion. The scheme as- 
sumes each material resides in its own sub-volume. We solve for separate, 
frequency-averaged radiation energy densities and matter temperatures in the 
sub-volumes. The energy change of the materials depends on the individual, 
frequency-averaged opacities as well as on AS. 

We now briefly describe the spatial discretization. We largely follow proce- 
dures described by Howell and Greenough [7] (H&G) and Shestakov et al [20]. 
Our MGD solver is embedded in an Eulcrian radiation-hydrodynamic code with 
cell-centered fundamental variables: p, u g , etc. The code has distinct 1, 2, and 
3D executables; mesh cells are line intervals, rectangles, and rectangular hexa- 
hedra, respectively. 
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In 2 and 3D, we discretize the diffusion term V • D g Vu g using the H&G 
subroutines since those are readily available. We note in passing that H&G use 
the Levermore-Pomraning flux limiter [10] instead of the simple expression in 
(3). For ID we have our own discretization; V • D g \7u g is written as 

[A+l/2(«i+l ~Ui)/h- A-l/2(Wj-U i _i)//l]//l, (15) 

where the group index g is suppressed and where i is the cell index. 

The /ace-centered diffusion coefficient D i+ i/ 2 is computed as follows. First, 
we modify (3) by adding the term (3/h to the denominator, where (3 is a small, 
user-specified constant, e.g., 10~ 6 . After factoring h, we obtain 

D = ch/[3 X h+\A(u)\/u + P], 

where we suppress the group index and note that the expression is to be evalu- 
ated on a face. The denominator is now dimensionless. The second term is the 
relative difference of u; we discuss its discretization momentarily. The product 
xh is an optical depth. In this light, (3 provides a floor to the cell's optical 
depth. The aim is to avoid complications with the matrix solve in case \ ~~ * 
and at the same time, u v is nearly spatially constant, which may easily happen 
for high frequency groups. The face-centered opacity is an average of the ad- 
joining cell-centered opacities. We offer several options. For the simulations in 
this paper, we typically use inverse averaging, but other options (arithmetric, 
square root) are also allowed. 3 The term \ A(u)\/u is written as 

2 \u i+ i - Ui\/(u i+1 + Ui) . 

Other options are also available, e.g., instead of the arithmetic average, one 
may substitute max(iii + i, Ui) in the denominator. We plan to extend the above 
discretization to higher dimensions. 

Cell-centered data, such clS Cy , SX6 obtained as in [7]. 

For coupling to the radiation field in mixed-material cells, we need averaged 
material properties, e.g., opacities. These are obtained by mass averaging. Sup- 
pressing the group index, if n is the material index and denoting averaging with 
an overbar, the opacity (cm 2 /g) is given by mR — ^2 n m n K n , where m n is the 
mass of the n th material. Equivalcntly, 

pK — ^ ' fnPn K n i 
n 

3 If the two opacities are very different, inverse averaging: K i+1/ / 2 = I {^i + — » 

2min(K i , K i+1 ). Assuming the opacity is monotone with T, the result is nearly the same as 
what is commonly done in gray diffusion, viz., forming a face-centered temperature, T i+ i/2 = 
(Ti +T;+i)/2, and calculating directly with T i+l / 2 - For example, if k = kq/T™ and 

Ti 2> inverse averaging gives 2ko/T™ while the face-centered T result is 2 n fco/T". 

For the free-free gray opacity, n = 3.5; hence, the two results are similar. Of course, if the 
opacity is not monotone with T, the face-centered technique is better. We plan to incorporate 
that option in the future. However, we note that multigroup opacities are usually not strong 
functions of T. 
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where /„ = V n /V is the volume fraction. 

This concludes the description of the algorithm used to advance the multi- 
group equations on an AMR level. In the following section, we analyze the 
convergence of (12)-(13), and we focus on how the ^tc parameter a ensures 
stable, robust iterations, to yield a physical, i.e., nonnegative result. 



3 Analysis of V&tc 

In this section we develop three criteria that set the \I/tc parameter. Disinter- 
ested readers can safely skip the analysis and continue to section 4 where we 
discuss the implementation of the multigroup scheme in the context of AMR. 

Recalling that a = 1 + t, we develop lemmas that set the initial magnitude 
of r, where by initial we mean the following. A new value of r is determined 
at each time advance for each AMR level. The level advance consists of nested 
loops. For the "inner" iterations, r is fixed. After convergence, r is reset to 
r — > a T r, where a T is a user-set input whose default value is 1/2. Section 5.2 
describes an experiment with another setting of a T . Our strategy for the initial 
t is to ensure a nonnegative w, diagonal dominance, and a convergent inner 
iteration. For the derivation, it is convenient to define 

G G 

B = Y,atB t , B' = Y,^B[, (16) 
C' g = a g B' g /pc v , C = B/pc v , C' = B'/pc v . (17) 



3.1 Positivity of w 

Before analyzing the effect of iPtc, we examine the scheme's behavior without 
it. If a = 1, the term u* disappears from (9). In the following discussion, 
we ignore the T° — T* term since for the SI scheme, or for the first FI inner 
iteration, T* =T°. Since vP g ~ B g , if cither At is large or the coupling is strong, 
a g B g ^> Ug. Hence, in this case, the RS of the system, w w a g B g — f g B, where 
Bis defined in (16). If a = 1, f g = a g B' g /(pc v + B') = C' g /(1 + C). Hence, 

w^a g (B g + B g C'-B' g C)/(l + C) . 

Since C and C are proportional to At times the opacity, the sole B g term in the 
numerator is swamped by the other two terms when At is large or the matter 
is optically thick. In this limit, the sign of w equals the sign of (B g C — B' g C), 
which may be negative. 

However, with \I/tc, nonnegativity of w is equivalent to the inequality 

< p(o-) = u* g a 2 + 2ba + c , 

where 

26 = u° g - u* g + a g B g + C u* g , 
c = C' (v,g — u* + a g B g ) + a g B' g [T° — T* — C] , 
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for the fully- implicit (FI) scheme. The SI scheme, adds the term a g B' (T* — T°) 
to the definition of 2 6. If a = 1, we recover the non-^tc scheme, which as shown, 
may have indeterminate sign(u>). At the end of the section we show that the SI 
scheme is less robust. We first analyze the FI scheme. 

For large a, p is positive — even if u* = 0. The derivative dp/do- increases 
monotonically and is positive for a = 1. If u* = 0, p increases linearly with a 
and has slope u° g + a g B g > 0. Hence, we have proved: 

Lemma 1 If p\ a =i > 0, the RS of (9) is nonnegative for all a > 1. Otherwise, 
(1) If Uj > 0, the RS of (9) is nonnegative if 



0~ > CTmin = HiaX 



(2) If u* — 0, the RS of (9) is nonnegative if a > cr mirl = — max(c/2 6). | 

The lemma's limit is very restrictive for large At, as we now show. As 
At — > oo, the terms a g , C and C dominate the definitions of b and c. Hence, 

lim 2 b = a„ B„ + C u* , 
lim c = a a B a C' — a a B'C . 

Substituting into the expression for the root and factoring out a g u* yields 

lim cr min = max^ (yj{a - (3) 2 + 4 7 - (a + (3)) , 

At— >oo Z \ / 

where a = B g /u*, (3 = Y,e K e,g B e/P c v, 7 = (B' g /u g ) J2e K i,gB' e / pc v and K£ i9 = 
Kg/Kg. The term a g — cAt ptt g equals £ c /£ g , where l c is the maximum distance 
a photon can travel in time At and l g is the absorption mean free path for 
the 5th group. We now show the remaining expression is of order one. If the 
radiation field is at equilibrium, a = 1. The term B' e is of order Bg/T. If it is 
exactly equal to Bg/T, the expression multiplying a g /2 vanishes. 

If u* = u° g = and p| CT= i < 0, then for large At, cr min — > (B' g C - B g C')/B g , 
which equals cAt times a term of order one. 

We now consider the SI scheme. As noted above, SI adds the expression 
a g B' g (T* — T°) to the definition of 2 b. Effectively, the extra term means that 
rather than having 2 b depend on the emission source B g (which is evaluated 
at T*), the coefficient depends on the linearization B g + B' g (T* — T ), with B g 
and B' g evaluated at T°. If the temperature is decreasing the expression may be 
negative. As a consequence, we are not assured that dp/ da is positive. If u* is 
nonzero, we can find a suitable a. However, if u* = 0, p(cr) is a linear function 
with possibly a negative derivative. If that case arises as we query the cells, we 
set a = 1 for the cell in question. Because of these uncertainties, by default, we 
run with the FI scheme. 
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3.2 Diagonal dominance 

To prove diagonal dominance, we compute row sums. The diffusion terms sum 
to zero, since the matrix composed of just these terms must annihilate the vector 
(1, 1, . . .). 4 Thus, for diagonal dominance, 

G 

a + a g - fg^di > d > . 
e=i 

Recalling the definition of f g , the relation is equivalent to 

< q{<y) = a 2 + 2ba + c , 

where 

2b = a g +C'-d, 
G 

c = a g C -C'g^at-C d, 
l=\ 

and C, C' g are defined in (17). As before, a > 1 is the domain of interest. The 
quadratic q(a) is nonnegative for sufficiently large a. However, 

G 

q\a=i = (l+C')(l + a g -d)-C' g 

The a g and C terms are proportional to At. Hence, as At — > oo, the sign of the 
expression is dominated by sign(C — C' g Y^eLi a e)- Since the expression varies as 
J2f=i a e(B'i — B' g ), the sign is indeterminate. However, (dq/dcr)\ a =\ is positive 
for d < 2. We have proved: 

Lemma 2 If q\ a=1 > and d > 0, A is strictly diagonally dominant for all 
u > 1. Otherwise, A is strictly diagonally dominant if 



O- > (Tmin = V b 2 - C - b . | 

Remark For large At, 

lim ( 7 min = max^ U {I - (3) 2 + A5 - {I + P)) , 

At— >oo Z \ / 

where S = {B' g /pc v )^2 e Ke tg and, as before, f3 = K£ ig B' e / pc v , and K£ g = 
Ke/ng. As in Lemma 1, when At is large, cr m i n = £ c /i g times a term which 
should be of order one. 

4 In extreme cases, because of finite precision, the diffusion terms may swamp the other 
terms. We discuss the possibility in Appendix B. 
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3.3 Two-step iterative scheme 

We've shown that for sufficiently large a, A is an M-matrix. Hence, (A — Mi) — 
M 2 and (A — M 2 ) — Mi are regular splittings, and each half of the two-step 
scheme (12)-(13) is a convergent iteration [25], Thm. 3.13, p. 89. Here we 
analyze how the scheme reduces the error. Of particular interest is that for 
large At, the scheme (12)-(13) may not converge unless the iPtc parameter a is 
sufficiently large. 

It is convenient to change variables, 



The system of interest is then A'v = w, where 

A' = A - M 1 - M 2 

and A is diagonal, 

A 9 = {T> d>g /a g ) -f g + l + oja g 

(Ml V)g = V .g Vg/a g 

G 

(M 2 v) g = f g J2 v e- 

If = v — defines the error for (12)— (13) , the error satisfies 

(A-Mi)e( J+1 / 2 ) = M 2 eW 
(A-M 2 )e( l+1 ) = M 1 e { - i+1 ' 2S > . (18) 

We express the error as a product of spatial and frequency components. For a 
2D spatial domain, 

4^ = ^e V=T(M *- H " 9 " ) . (19) 
where the indices k and m refer to distinct spatial axes. We now analyze the 
iteration error 

e(* +1 > = (A - M 2 y x Mi (A - Mi) -1 M 2 e w . 

Consider the unit vector consisting of N components, with unity in the £th 
position and zeros for the rest. Since the initial error is a linear combination 
of such vectors, it suffices to analyze the case when the frequency component of 
eW equals e>. We will prove that for a properly chosen a > 1, 

||e< i+1 >||i<C<l. 

In other words, if a is sufficiently large, one iteration of the two-step scheme 
reduces the error, which we will show occurs for (' < (, cr(C') > c(C)- However, 
larger a denote a smaller \Ptc time step, resulting in a longer pseudo-time to 
reach the desired steady-state. 
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Assuming that the diffusion coefficient does not vary in space 5 and that we 
use a uniform 2D spatial mesh with mesh size h, the error after the first half 
step is 

4 i+1/2) = fg 1 1 1 - fg + WS) + 2 Vg (2 - cos 9 k - cos m ) } , (20) 

where 

Vg = Dg/dgfl 2 = lg tD I g , a / 3h 2 , 

Ig.a — ^/P K g is the absorption mean free path, and l g ^u is the diffusion mean 
free path; the latter is the sum of the absorption and scattering opacities. In 
(20), the expression multiplying r\ g is nonnegative. 6 Since f g < 1, e g l+1 ^ 2 ^ is 
nonnegative. Assuming the worst case 9 — yields, 

0<e g i+1 ^<fg/(l-f g + a/a g ), (21) 

a result which also holds in 1 and 3 dimensions. The bound is sharp; i.e., e^ 1 / 2 ') 
equals the bound if the original error has no spatially varying component. 

Since (19) holds for i, i + 1/2, and i + 1, we now analyze the second half 
step. In 2 dimensions, 

Mi e (i+i/2) = (cog Qk + cog m )2r) g e (i+1/2) . 

In n = 1, 2 or 3 dimensions, the parenthetical expression contains 1, 2 or 3 
cosine terms. If we again assume 9 — 0, the expression is bounded by n. 

To determine from (18) we invert (A — M2) using the Sherman- 

Morrison formula by noting that 

A-M 2 = A'-fe T , 

where e is the vector consisting of all ones, the components of f are the previously 
defined f g , and A' is diagonal with 

M g = rf g + \ + a/a a , ^2n % . (22) 

In (22), we generalized by allowing for n = 1, 2, or 3 spatial dimensions. After 
some algebra, we obtain 



| c (i+i)| < — 

1 3 1 - A' 

9 



G , U+l/2) 



, (i+1/2) , f fg \ Vt^t 

1,1 •■' \ 1 - n '' [ \M f J ^ 



where 



l-e T (A')- 1 f =1-^/^. (23) 



5 Since, as wc show, the worst error arises for spatially constant error, we are free to ignore 
the diffusion flux limiter in the analysis. 

6 The corresponding expression in 1 and 3 spatial dimensions is also nonnegative and 
bounded by f .0 and 3.0, respectively. 
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Summing yields the 1-norm, 



G i (i+l/2) 

< [l-e T (A')- 1 f]" 1 Y^^ 1 
1=1 



A' 



Our task is done if we can show that the RS is bounded by (. Using (23) this 
entails showing that 

G „/, (i+l/2) 

E 

£=1 

After substituting the bound (21) and simplifying, the inequality becomes 




G 

E 

9=1 



a - 



C' 



G 



g=l V 9 



(24) 



To summarize, if (24) is satisfied the two-step scheme (12)-(13) converges and 
each iteration reduces the error by a factor (. 

We now show that if \ttc is not used, i.e., if a = 1, and At is large, the 
scheme may not converge. If a = 1, since a g oc At, liniAt^oo = rj' g + 1 and 
liniAt^oo f g = where p s > and X)oPs = 1- Also, if cr = 1, since C' g oc At, 
for large At, the lone a on the RS of (24) is swamped by the sum. Dividing 
both sides of (24) by cAt, the LS becomes 



E 

9=1 



pn g B' G r)' g 
pc v {l + r}') 



(1 -Pa)" 1 - 



On the other hand, if £ = 1, the RS tends to the same sum, but without the 
term (1 — Pg)^ 1 ■ This makes the LS larger than the RS, giving the desired 
contradiction. We have proved: 

Lemma 3 If a = 1 and At is large, (12)- (13) may not converge. \ 

We now estimate how large to make a in order to satisfy (24). The terms a g 
and C' g are proportional to At; also, a > 1 and A^ > 1. Hence, (24) holds for 
small At. To obtain a tractable expression, we derive a relation that stems from 
a more stringent inequality. Equation (24) holds if we derive a a that satisfies 
a relation insensitive to the lone a on the RS and is obtained by requiring that 
the individual terms in the sum satisfy the inequality. This allows canceling the 
common term C' g /A' g . Hence, we seek a satisfying 

rf g /(l-f g + a/a g )<C(A' g -l). 
Recalling that f g — C' g /(a + C) and using (22) leads to 

< s(a) = (T 3 + a s a 2 + (3 S a + 7s , 
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where 

a s = a g (l + rf g )+C, 

ft = a g [{l + V ' g )C' -q + agr/^l-C 1 )], 

As before, a > 1 is the domain of interest. 

To simplify the analysis, we assume ( = 1, i.e., we seek a a that guarantees 
marginal convergence. To this end, we define 

P s = a g [(l + rj' g )C'-C' g }, 
Is = -a 2 g V g C' g . 

Consider the cubic 



s(a) = a 3 + a s a 2 + f3 s a + j s . 
For a > 1, all derivatives of s are positive. If 

s(l) = l + a s + /3 s + 7 s >0, 

then the scheme (12)-(13) converges. However, if s(l) < 0, we need a a > 1 
that renders s > 0. To avoid computing cubic roots, we approximate s by a 
quadratic w(a), 

w{a) = (3 + a s ) a 1 + S - 3) a + 7, + 1 , 

and determine the root of w. The polynomials w and s and their first two 
derivatives agree at u = 1. The difference s(<r) — w(cr) = (cr — l) 3 , i.e., w(cr) < 
s(a) for ct > 1. Hence, the positive root of w(a) overestimates the a needed for 
marginal stability. We have proved: 

Lemma 4 If w\ a =i > 0, the scheme (12)-(13) converges. If w\ a= i < 0, the 
scheme converges if 

J 0s - 3) 2 - 4 (3 + a s ) (% + 1) + 3 - 0„ 

l + T = a> <r min = — — . I 

6 + 2a s I 



4 Multigroup AMR scheme 

In this section, we describe our implementation of AMR for the multigroup dif- 
fusion (MGD) system. The scheme necessarily adheres to the code's general 
architecture. That is, on each grid level each physics module (hydrodynam- 
ics, radiation) is called in order. These comprise the level solves. If AMR is 
used, the code refines in both space and time, as described by Howell and Gree- 
nough [7]. After a refined level is time-advanced to the next coarse level time, a 
synchronization is required in order to maintain conservation. For a scalar dif- 
fusion equation and only two levels, coarse and fine, the "sync-solve" is difficult 
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enough since it reduces to effectively an unstructured grid solve over the com- 
bined coarse and fine grids. For MGD, the difficulty is compounded by having 
to sync-solve a coupled system of diffusion equations. 
We begin by recalling the equations, 

d t Ug = V • DgX/Ug + Kg(Bg- U g), Q = 1, . . . , G (25) 

G 

c v d t T = -^KgiBg-Ug), (26) 

9=1 

where c v is now the heat capacity, while D g and the diffusion and cou- 

pling coefficients. For ease of exposition, it is convenient to consider the one- 
dimensional case. The level solve module computes the solution to 

« fl ,i-<i = {Fga+i/2-Fg^ 1 , 2 )/h l + lg , l [Bg{T l )~u g , l ], (27) 

G 

Cy^Tt-Tf) = -$> fl ,i[B fl (2i)-u fl ,i], (28) 

9=1 

where i is the cell index, 7^ = Atn 9i i, and F g i+1 / 2 is the fluence on the right 
edge of the ith cell, 

^5,1+1/2 = ^tD g i+1 / 2 (%,j+l - U g s)/hi . 

For simplicity, assume there are only two levels, coarse and fine. Since (25)- 
(26) are reaction-diffusion equations, advanced with backward Euler temporal 
differencing, the discretization is unconditionally stable. Hence, in the following, 
in order to simplify the derivation, we assume that both levels are advanced with 
the same timestep. However, in the code we also time-cycle. If i = 1, . . . , N 
define the indices of all coarse-level cells, let j = 1, . . . , J define the indices of 
the refined cells and i = I, . . . , N define the indices of those coarse cells which 
are not refined. Coarse cells indexed with i — 1, . . . ,1 — 1 arc defined as the 
"covered" cells. We first update the entire coarse level, then the fine level. 
Both levels require boundary conditions (BC). The coarse level uses the user- 
specified BC. In the following example, the refined domain abuts the left side 
boundary and consists of J cells. Hence, the fine level uses the same BC on the 
left edge. The fine cell indexed with j = J lies in the interior of the domain. 
We reuse the Howell and Greenough [7] infrastructure to provide a Dirichlct 
condition for the cell. The datum is obtained by interpolating coarse grid data. 
Let kj and hi define the mesh widths of the fine and coarse cells, respectively. 
After multiplying by the mesh widths and summing over all cells and groups, 
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we obtain 



3 = 1 



g=i 



N 



i=J 



9=1 



E ( F 9,N+\/2 ~ ^g,l/2 - ^,c/) , 
9=1 



(29) 



where the last term is the fluence miss-match of the gth group at the coarse-fine 
interface, 



SF, 



g,cf 



g, .7+1/2 



F g,i-i/2 



The AMR scheme assumes that the system is linear. Hence, the emission is 
expressed as 

B g (T i ) = B- g . + B;.(T i -T-), 

where B* i and its derivative with respect to T, i.e., B gi , are evaluated at 
a temperature T*, e.g., B* i = B g (T'). For semi-implicit Eulcr differencing, 
T* = Tf; if fully implicit, T* = Ti. Either way, because 5F g ^ c f need not be 
zero, (29) shows that energy may not be conserved after the two level advances. 
To restore conservation, we introduce the system for the corrections 



9,« 



c T' 



( F Li+i/2 ~ F 'g,i-i/2)/hi + 79,, [B' g ,i Ti - u'J + 5F g>cJ /hi , (30) 



G 

E- 

9=1 



[ B g,i T i 



(31) 



where F' g i±l/2 



denote the implicit fluxes; they are functions of u' g . 
Equation (30) holds for all groups g = 1, . . . , G. In (30)-(31), the mesh 
index i varies over the coarse cells not marked for refinement (i = I,... ,N) as 
well as the fine cells (j — 1 , . . . , J) . Following the methodology of [7] , we put 
the fluence mis-match SF g ^ c f into the coarse cell(s) abutting the interface of the 
coarse and fine domains. 

Summing the level-advance and correction solutions yields conservation. If 
u * g ,t = u g,i + u ' g .i and T* =% + T-, combining (30)-(31) with (27)-(28), mul- 
tiplying by the mesh widths, and summing over cells and groups, yields the 
desired conservation relation, 



3=1 



9=1 



+ 



N 



i=i 



Cv ,(T;-Tf) + J2Ki-<i) 

9=1 



= E { F g,N+l/2 _ F *g.i/2) ■ 



9=1 
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Equations (30)-(31) present a formidable task as it requires solving a si- 
multaneous system of equations for (G + l)N unknowns, where N denotes the 
number of refined cells plus the number of coarse cells not marked for refine- 
ment. The grid is effectively unstructured since it combines coarse and fine 
discretizations of the domain. We attack the problem by applying a variant 
of the "Partial Temperature" scheme [11], [19]. In this scheme, groups are as- 
signed a random order. As we cycle through the groups, each group computes 
a correction u' g and a partial temperature T g . Note the group index g for the 
temperature. Although the scheme decouples the groups from each other, the 
partial temperature T g changes as we cycle through the groups. To be precise, 
for each group, we solve the system 

U 'g,i = ( F 'g,i+l/2 - F 'g,i-l/l)/hi + 

JaAB'g^i-u'J+SF^f/hi, (32) 
Cv,i (T g ,i - T fl _ M ) = -j 9yl [B' t T g4 - u' g i ] , (33) 

where, as above, the mesh index i ranges over all refined cells and all coarse 
cells not covered by the fine grid. For the group index g\ that we first pick, 
T g -i,i — on the LS of (33). Solving (32)-(33) for g = g\ yields the first partial 
temperature T gi . This temperature replaces T g _i.j on the LS of (33) for the 
second randomly picked group gi. After cycling through the groups, the last 
one, go, gives the desired corrected temperature, i.e., T[ = T' ga i . 
If (33) is summed over all g, the LS telescopes and we obtain, 

G 

Cv,i Ti = — "fg^i [B g i T g j — u g i ] . 

3=1 

Because we have T Qyi on the RS instead of T/, this is not exactly (31). However, if 
T g ,i doesn't vary too much as we cycle through the groups, the result is no worse 
than one obtained with the (commonly-used) partial temperature (PT) scheme 
since we apply PT to only corrections of the level-solve solution. Cycling through 
the groups in random order avoids biasing the deviation since the coupling in 
(32)-(33) may lower T for one group while raising it for another. In any case, 
the combined solution (u*, T*) is still conservative. 

Equations (32)-(33) are solved using a Schur complement. Since (33) does 
not involve spatial derivatives, we can easily solve for T^. After substituting 
the result into (32), we obtain a single scalar equation for u' g i , albeit now, on 
the unstructured grid composed of coarse and fine cells, viz., 

= 0^+1/2 - F 'g,i-l/2)/hi + 7 9 ,* Vg,i T 9~hi - u 'gJ + ^ F 9,cf/hi , 

where, r/ Si j = c v ^ / (c Vyi + 1 g ,i B g i ). After solving for u' g i , equation (32) yields 
T g ,i. The fluence miss-match SF 9tC f acts as a source to the corrections. For 
groups with long mean free paths (mfp) and weak coupling, 5F g _ c f diffuses over 
the mesh. For groups with short mfp and strong coupling, 6F g , c f is spread 
locally over the group energy u' g t and "absorbed" into the matter. 
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Before closing this section, we note an inconsistency in the above multilevel 
scheme, indeed in any scheme embedded in a multi-physics code like ours, which 
advances several modules (hydrodynamics, heat conduction, radiation) using 
operator splitting. With splitting, on each level, the modules are advanced in 
order. For simulations using hydrodynamics and radiation diffusion and running 
with coarse LO and fine LI levels, the order of operations is as follows. Level 
LO first advances hydrodynamics, then radiation. Next, if refining by a factor 
of two, LI advances in the order: hydrodynamics, radiation, hydrodynamics, 
radiation. The multilevel solve advances in the same order: hydrodynamics, 
then radiation. This implies that the radiation multilevel solve uses coefficients, 
e.g., p, that arc not the same as those used by the radiation level solve mod- 
ules. In principle, one cannot simply add the correction equations to the level 
solve equations and claim that the sum satisfies a consistent set of equations. 
Nonetheless, the solution remains conservative. 

5 Simulations 

This section presents results using the multigroup scheme. We consider three 
problems. In Section 5.1, we present a test problem with a known analytic 
solution. We compare numerical results with tabular data, previously published 
by Shestakov and Bolstad [22]. Using Richardson extrapolation, we show that 
our \I>tc scheme, i.e., what we apply on a level, is second (first) order correct in 
space (time). When running with AMR, the temporal accuracy is first order. 
Accuracy of the spatial order depends on the norm used to measure convergence. 
In the most stringent oo-norm, the order degrades to first, or worse, as shown 
at the end of section 5.1. Section 5.2 develops a variation of the Section 5.1 
test problem in order to demonstrate the benefits brought by \I/tc. We do 
this by running with and without \ftc. We make several runs, each for only 
one timestep. Runs are made with successively larger At. Because of fully 
implicit differencing, as At — > oo, the numerical solution should approach the 
time-independent, steady-state. The problem in Section 5.3 brings everything 
together. We simulate the explosive expansion of a metal sphere suspended in 
air. The expansion is due to sourcing a large amount of energy in a short time 
into the sphere. Simulations are done with the code's full functionality, i.e., we 
couple all of the physics modules and also use AMR. 

5.1 Linear MGD test problem 

In this section we present results for a MGD problem with a known solution. 
Due to the nonlinearity of the equations, there are no test problems with an- 
alytic solutions. Thus, to validate and verify our algorithm, we consider the 
linearized multigroup equations developed by Shestakov and Bolstad (S&B) 
[22] and compare with tabular data. 

The S&B tables present results for a 64-group discretization of the linearized, 
nondimcnsional, multifrcquency diffusion equations derived by Hald and Shes- 
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takov (H&S) [6]. In the following, we briefly derive the nondimcnsional system, 
describe the test problem, explain how to set up the problem in a radiation- 
hydrodynamic code, demonstrate the problem's relevance to typical applica- 
tions of multigroup diffusion, compare results with an improved-accuracy table 
[4] (supplied in the Appendix), and conclude by proving that our multigroup 
scheme's convergence is first order in time and second order in space. 

The nonlinear multifrcquency H&S system is derived by assuming slab sym- 
metry, constant density, an ideal gas EOS, and an opacity characteristic of 
free-free transitions. One advantage of the H&S system is its nondimensional 
form, which enables comparing results from codes using different dimensional 
units. The equations are obtained by choosing characteristic values for den- 
sity po, temperature T , and inverse mean free path (mfp) k = k /v 3 with 
kq = const and v the frequency variable. Radiation emission is given by a Wicn 
distribution 7 , i.e., Bw = B v 3 exp(— hv/kT), where B = 8n h/c 3 is the same 
constant defining the Planck function. The inverse mfp appears in both the dif- 
fusion, D = c/3k, and the radiation-matter coupling terms, en. (The diffusion 
is not flux- limited.) The normalization proceeds as follows. The values po, kq, 
and T define the other normalization constants, 

is = kT /h, £ = v^/k q , x =£ /\/3, 
h =£o/c, u = B vl , E a = u v . 

By defining nondimcnsional variables, x' — x/xq, t' = t/to, u 1 = u/uq, v' = 
v/vo, etc., (and dropping the primes) we obtain the normalized system, 8 

d t u = V -v 3 Vu+ {v 3 e- v ' T -u) / v 3 , (34) 

p CO 

Rd t T = -T+ (u/v 3 )dv, (35) 
Jo 

where the constant 

R = (h/k) (poc v /u a ) 

and c v is the specific heat. Henceforth, unless stated otherwise, we use nondi- 
mcnsional variables. 

The H&S system yields a precise definition of the multigroup equations since 
the group integrals can be computed exactly, an impossible task for definite inte- 
grals of the Planck function. Given a group structure {v g }f=o, after integrating 

7 It is noteworthy that H&S's choice of opacity and Wien spectrum for B gives the same 
emission source n Bw as would be obtained by including stimulated emission (SE) effects 
[26] and using the Planck function, since SE multiplies k by the factor (1 — e~ hv / kT ). Also 
note that without SE, the resulting Planck-averaged gray opacity does not exist; the integral 
diverges. 

8 If instead of Byy, H&S had used the Planck function, the factor e~ v l T in Eq. (34) would 
be replaced by (e v / T — l) -1 . However, H&S would then be unable to form Eq. (35), since 
the integral over all v (the total emission) diverges — see prior footnote. 
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over groups, 

d t u g = v 3 g d xx u g +p g T -u g /v 3 g , g=l,...,G (36) 

G 

Rd t T = -T+Y^ug/P* (37) 

9=1 

where u g — J g u dv and v g is a group's representative frequency. S&B define 
v g as ^]VgVg-\ and V\ as v\/2 since the lowest group boundary is zero. The 
emission coefficients are 

p g = exp(-f s _i/T) - cxp(-^ g /T) . (38) 

If the group structure is broad enough, ^2 g p g = 1. 

Equations (36)-(37) are nonlinear because of the product p g T. To derive 
an analytic solution, S&B follow the approach of Su and Olson [23], [24], which 
requires a linear system since it uses Fourier and Laplace transforms. S&B 
linearize by defining a fixed temperature Tf and substituting Tf for T in (38). 

Except for one item, it is easy to assemble the S&B linearized MGD sys- 
tem in a conventional radiation-hydrodynamic code. Such codes usually allow 
an ideal gas EOS and a desired analytic form for the opacity. One chooses 
arbitrary values for po, no, Tq, and picks a specific heat c v to set R. In our sim- 
ulations, pa = 1.8212111 • lCT 5 g cm" 3 , T = 0.1 keV, and n = 4.0628337 • 10 43 
cm -1 s -3 . To comply with S&B, we chose c v to obtain R=l. Our po, T , and 
kq choices were dictated purely by reasons of convenience. Since we compare 
with a nondimcnsional result, other constants may be used instead. 

The subtle item is how to force a code's spectral emission rate to equal 
p g (Tf)T. We accomplish the task as follows. The gth group's emission is 
a g [Bg+B' g ( T—T* )], where a g = At cp n g and n g is the group-averaged opacity. 
The terms B g and B' are integrals over the gth group, at temperature T* , of 
the Planck function and its derivative w.r.t. T. The integrals are computed by 
a FORTRAN subroutine, which takes T* as an input variable. For the test 
problem, we use a different subroutine, which when called, first defines 

B 'g = (^o) 3 (87rfc/c 3 )[exp(-y 3 _ 1 ) - exp(-y fl )] , 

where y g = hv g VQjkTjT§. After computing B' , the routine sets B g = B' g T* . 
In the y g definition, v g and Tf are nondimcnsional, while v and T are the 
normalization constants. The (DgVo) 3 term cancels the \jv z dependence of the 
opacity. 

For the test, we consider S&B's problem 1. The nondimcnsional domain 
is < x < X, where we set X = 4. The initial condition is T = 1(0) for 
x < (>) 0.5 and u = everywhere. We use symmetry boundary conditions at 
x = and homogeneous Milne at x = X, i.e., u g + (2£ g /3) d x u g = 0, where £ g 
is the mean free path. We use the same group structure as S&B: 64 groups, 
starting at zero, with widths increasing geometrically by the factor 1.1. We set 
v\ = 5 • 10~ 4 as the width of the first group. 9 The test simulates an initially hot 

misprint in [22] erroneously has v\ = 10 . 
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Figure 1: Linear MGD test. Comparison of the linear solution (Tf — 1.0) with 
the solution of the nonlinear MGD system with Planckian emission; t — 1. 



slab of material encased by cold matter. Since u is initially zero throughout, the 
solution evolves by first coupling in the hot subdomain. As radiation diffuses 
out, it couples to cold matter thereby heating it. Because of the opacity's l/i/ 3 
dependence, the group's diffusion and coupling rates differ. 

Although the problem appears contrived, it represents effects of radiation 
diffusion. We prove the assertion in Fig. 1 where we display the temperature T 
and the total radiation energy density E r (= u g ) for two simulations ending 
at t = 1. Solid lines pertain to the linearized system, where Tt = 1.0. Dashed 
lines are solutions of the "physical" nonlinear MGD system using Planckian 
emission. The similarity of the solutions validates the relevance of the test 
problem. We used Tf — 1.0 (instead of S&B's Tt — 0.1) because over the short 
duration of the simulation, the emission temperature in the hot subdomain is 
of order 1.0 rather than 0.1. 



5 SIMULATIONS 25 



X 


e(T) ■ 10 a 


s(E r ) ■ 10 3 
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e(T) ■ 10 a 


e(E r ) ■ 10 3 


0.00 


0.0016 


0.3012 


0.51 


4.8468 


0.2785 


0.20 


0.0015 


0.3028 


0.52 


1.8220 


0.0031 


0.40 


0.0005 


0.3268 


0.53 


1.0528 


0.1293 


0.46 


0.0081 


0.3903 


0.54 


0.7320 


0.2128 


0.47 


0.0174 


0.4252 


0.60 


0.3316 


0.5263 


0.48 


0.0467 


0.4945 


0.80 


0.6099 


1.3841 


0.49 


0.2205 


0.6979 


1.00 


1.4253 


2.2138 


0.50 


0.0019 


0.3518 









Table 1: Linear MGD test. Relative errors times 1000. Numerical result ob- 
tained with T f = 0.1, h = 1/400, At = 1/200. 



We now present our MGD result using S&B's parameter Tf — 0.1. Table 1 
displays the relative errors of T and E r for various x, at t = 1.0. For a variable 
/, we define the error e(f) = \(f x — fk)/fx\, where fk are our numerical results 
and f x are the S&B table values, listed in the Appendix. 10 Table 1 shows 
that we obtain better than 0.5% accuracy over the domain < x < 1. The 
worst error 0.48% occurs for T at x = 0.51. At that point, according to the 
table in the Appendix, T undergoes more than a 20-fold drop from its value at 
x = 0.49. We focus attention at the domain near x — 0.5 since that is where the 
variables undergo the sharpest change. At these points, we obtain better than 
0.1% errors, except for T at x = 0.52 and 0.53. Errors near x = 1.0 are less 
important for two reasons. First, the S&B domain extends to infinity while ours 
extends to only X = 4. Hence, for large x, our results become less accurate. 11 
Second, our code requires having a positive min(T). Hence, we cannot initialize 
with T — in the cold region. At the end of the run, at x = 1, our temperature 
has risen by only a factor of 10 4 , which precludes reaching much better than 
0.1% accuracy there. 

We were unable to use the S&B tables for a convergence study to verify 
our scheme's convergence properties w.r.t. timestep At and mesh size h. We 
speculate that the reason is that the truncation is a mix of errors due to finite 
At and h. Hence, a refinement study of one may be polluted by an overly coarse 
value for the other. However, we can use Richardson extrapolation to prove that 
our scheme is correct to first order in time and second order in space. Let Vk 
denote a numerical solution to an equation discrctized by a constant parameter 
k. For an initial value ODE, k represents the timestep; for a time independent 

10 For each point, /j. is the arithmetic average of the two adjoining cell-centered values. 

n When we compare results of two simulations at the cells adjoining x = 1.0 where one run 
uses X = 4 and for the other, X = 8, wc find the relative differences: 8 ■ 10~ 6 , 2 • 10~ 7 for 
E r , T, respectively. Since these differences are 3-4 orders of magnitude less than the Table 1 
errors at x = 1.0, increasing the domain beyond X = 4 would have little impact on the entries 
of Table 1. 
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Figure 2: Timestep and meshsize orders of convergence; At (Ax) on left (right) 
sides; t = 1.0; see text. 



equation, k is the mesh width. If v is the analytic solution, 

v k = v + ak a + 0{k b ), 

where < a < b, and where a is independent of k. In the asymptotic regime, the 
k a term dominates the error, which allows ignoring the 0(k b ) term. Assuming 
we have three solutions Vk, W2fc, ^4fc, a ratio of differences yields 

v 2k - V4k = 2 a 
Vk - V 2 k 

The order of convergence a is found by taking logarithms. 

We apply this procedure to estimate the orders of convergence. First, for 
the At study, we fix h = 0.01 and obtain three results using k = 0.5 • 10~ 8 s, 
2k and 4k. For the Ax study, we fix At = 0.5 • 10~ 8 s and use k = 0.0025. In 
both studies, runs are halted when t — t - We compute a at 15 points across 
the domain [0, 1] for both E r and T and focus attention at x = 0.5, where the 
fields undergo the sharpest change. Results are presented in Fig. 2. The left plot 
clearly displays first order temporal convergence since a ~ 1 across the domain. 
The right plot supports our contention of second order spatial convergence. The 
low a ~ 1.82 (1.89) values for E r (T) arise only at the two points x = 0.49, 0.51. 
We claim that at these points, we are not yet in the asymptotic regime. 

The results of Fig. 2 pertain to a solution obtained on a single level, i.e., 
without using AMR. We now analyze how AMR affects the order of spatial and 
temporal convergence. For each study, Ax and At, we make three simulations 
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Figure 3: Timestep order of convergence on composite AMR mesh; see text. 



(as before, we halt at t = 1.0) in order to apply our Richardson extrapolation 
technique. In each study, the composite grid consists of a "base" level L0 mesh 
over the entire domain and two AMR levels. Each level refines by a factor of 
two. Both LI and L2 levels refine around x = 0.5. We examine convergence at 
points x in all levels. 

For the At study, all three runs use the same composite spatial mesh. We 
make three runs; each with fixed timesteps Ato, 2 Ato and 4Ai , where At = 
1/400. The composite mesh uses Ax = 1/100 on level L0. The LI mesh extends 
over 0.36 < x < 0.64, and the L2 mesh extends over 0.42 < x < 0.58. We obtain 
nearly the same temporal order as for the level solve. Figure 3 displays a for the 
128 cells on < x < 1. The lowest order, a « 0.94, occurs at x = 0.38 (0.62) 
for E r (T) near the L0 and LI coarse-fine interface. 

The Ax study requires more care. For each run, the LI mesh extends over 
0.25 < x < 0.75, and the L2 mesh extends over 0.375 < x < 0.625. We refer 
to the three runs as Rl, R2 and R4, where Rl and R4 use the "coarsest" and 
"finest" composite grids, respectively. For the three runs, the level L0 mesh 
sizes are 1/40, 1/80 and 1/160, respectively. Because each AMR level refines 
by a factor of two, for Rl, the L0, LI and L2 mesh sizes arc also 1/40, 1/80 
and 1/160. The R2 mesh widths are 1/80, 1/160, and 1/320; R4's are 1/160, 
1/320 and 1/640. The composite grids are constructed so that within each level, 
the Rl cell boundaries are also cell boundaries of runs R2 and R4. Hence, by 
arithmetic averaging adjoining cell-centered data, we obtain numerical results at 
the same points for each run. These (averaged) values are used for Richardson 
extrapolation. Figure 4 displays the ratio {fm — fRi)/{fm — Ir2) for the 79 
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Figure 4: Meshsize order of convergence on composite AMR mesh; see text. 



faces on < x < 1. The ratio is approximately 4 over most of the domain, which 
indicates second order convergence. However at the coarse-fine interfaces, the 
order drops significantly; especially for T at x — 0.25 and 0.375. 

The loss of accuracy at the coarse-fine (C-F) interfaces is due to the dis- 
cretization of the diffusion operator. We use the infrastructure developed by 
Howell and Greenough [7] to assemble the linear systems. Unfortunately, the 
difference stencils — which are not discussed in detail in [7] — have a shortcoming 
near the interface. A more accurate discretization would yield an asymmetric 
matrix; for reasons of efficiency, symmetric linear solvers were preferred. 

The inaccuracy can be analyzed by considering a derivative such as u xx near 
the C-F interface. Assume that level L0 lies to the left of LI. For i = 0, 1, 2, let 
Xi denote the first three cell centers on LI and let h define the LI mesh size. Let 
x c denote the center of the coarse cell next to the C-F interface. On LI interior 
points, e.g., on x 1; u xx is approximated by the difference: (u — 2u x + u 2 )/h 2 . 
Hence, 1/h 2 is the off-diagonal matrix coefficient corresponding to u on the 
x\ row. For the matrix to remain symmetric, the u\ coefficient on the xq row 
must equal 1/h 2 . At xq, u xx is written as a difference of the right and left 
fluxes divided by the cell width h. The right flux is (ui — u )/h. The left flux 
is expressed as the difference (u — u c ) divided by the distance between the cell 
centers. If LI refines by a factor of two, the distance x — x c = 3h/2. Thus, at 
x , to maintain symmetry, u xx is approximated by 



/ Ul - Up _ u - u c \ I 

\ h 3h/2 ) I 
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Unfortunately, the left flux is not centered on Li's left-most face (at x = xo — 
h/2). A Taylor expansion shows that the difference is inconsistent; it equals 
(5/4:)u xx + 0{h) and this is the source of the error. However, the error is 
localized. In a global sense, it is O(h), when computed by integrating over 
the entire domain: J u xx dV. This concludes the refinement study on an AMR 
mesh. 

To summarize, in this section we have shown: (a) With a proper choice of 
Tf, the test problem mimics MGD physics, (b) We obtain excellent agreement 
with the S&B tables, (c) Our scheme is correct to first order in time and second 
order in space, (d) On an AMR mesh, the scheme incurs the same loss of 
accuracy as the one presented by Howell and Grccnough [7] since we use the 
same discretization at coarse-fine interfaces. 

5.2 Benefits of *tc 

We now present results that illustrate the benefits obtained by using \I/tc. We 
show that for sufficiently large At, the conventional (ADR) scheme of Axelrod 
et al [2] i.e., where a = 1, fails to converge. Furthermore, if At is only moder- 
ately large, so that the ADR scheme does converge, introducing vl/tc accelerates 
convergence. 

We begin by considering a variation of the problem introduced in Section 5.1. 
In this section, unless stated otherwise, we use normalized variables. First, we 
replace the Wien distribution with the Planck function. After normalizing, we 
obtain an equation similar to (34) except that e~ v l T is replaced by {e u / T — 
Without stimulated emission effects, the multifrequency system is ill posed since 
the RS of the temperature equation integrates the coupling term over all v. (The 
integral of B v jv z diverges.) Since this is only a test, we ignore this complication. 
We use seven geometrically spaced groups, whose widths double with increasing 
frequency. The leftmost group boundary is zero; the first group width V\ — 0.5; 
the last boundary 1/7 — 63.5. As in Section 5.1, the first group's opacity is 
evaluated at v\j2 and the rest are evaluated at the square root average. The 
spatial domain is < x < 2. The initial conditions are as before, viz., T = 
1 (0) for x < (>) 0.5 and u is initially zero. We impose symmetry boundary 
conditions on both left and right endpoints. Hence, at all times, the total energy 

1 /2 

should equal the initial amount J Q RT dx = 1/2. 

Our test consists of several runs, each for only one timestep. All runs use h = 
0.01. We run in fully implicit mode; hence, upon convergence, the temperature 
T and emission source B V (T) are consistent. For infinitely large At, a single 
time advance yields the steady-state with T = T r , where the radiation energy 
E r = aT^. In the nondimensional system, since B v is the Planck function, 
a = 7r 4 /15. Hence, the equilibrium temperature is the solution to, 

2(T e + aT e 4 ) = 1/2, i.e., T e = 0.2314. 

The ^tc result, where At — 1000, is displayed in Fig. 5. The figure shows that 
the two fields are nearly in equilibrium and almost spatially constant; T r and 
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Figure 5: ^tc robustness test; solution after one time advance; h = 0.01, At = 
1000. 



T vary less than 1% and 2.4% respectively. The initially high T in x < 0.5 has 
decayed more than fourfold. The radiation field, as it coupled in the initially 
hot region, diffused outwards thereby heating the cold region. 

The simulations were run with and without \I/tc. Both runs consist of nested 
"inner" and "outer" loops. The inner iterations (12)-(13) progress until the 
residual and the iterate difference | |it( i+1 / 2 ) — 1 1 fall below specified tolerances 
(which may not happen). At that point, the outer iteration computes a revised 
temperature T using (6). We then reset T* = T and use it to recompute the 
B( and B' e coefficients. For the first outer iteration, T* — T°. The iterations 
conclude when the temperature change and the nonlinear residual fall below 
their specified tolerances. 

The problem's difficulty increases with At. Without ^tc, it becomes im- 
possible to solve if At is very large because of the computer's finite precision. 
For large At, the time derivatives, e.g., (it — u°)/At, are dominated by the 
other terms. Hence, the initial condition (u°, T°) becomes less relevant. Un- 
fortunately, energy conservation depends on "remembering" the initial condi- 
tion. The boundary conditions enhance the difficulty. If the initial condition 
is indeed "forgotten," the solution is not unique. Any equilibrium temperature 
T e = T = T r is a steady-state. 

For runs without \&tc, we impose a = 1 and determine for which magnitude 
At the iterations fail to converge. Runs using ^tc proceed as follows. We 
first compute the three different a required to have (1) a nonnegative RS, (2) 
diagonal dominance, and (3) convergence of inner iterations. That is, the a must 
satisfy the lemmas of Sections 3.1, 3.2 and 3.3. The iterations commence using 



5 SIMULATIONS 



31 



the largest a. The parameter a is fixed for each outer iteration. Experience 
has shown that the lemmas give an overly large a. Hence, we use the lemmas 
to set a for only the first outer iteration. Subsequent outer iterations decrease 
a as follows. Recall a = 1 + t and that only when r = do we solve the 
correct discretization of the equations. Successive outer iterations multiply r 
by a constant factor, i.e., r — > a T r. The factor may be changed by the user. 
For small a T , r decreases quickly, but the resulting linear system is harder to 
solve. For the hardest test, where At = 1000 (see below), we experimented and 
found better results with a T = 0.5 than with a T = 0.25. 

Our tests begin with At = 20, a magnitude at which both modes, with and 
without \I/tc, converge and give nearly identical results. For this moderately 
large At, \ttc brings the benefit of faster convergence: 37 vs. 50 CPU sec, i.e., 
nearly 33% faster. For At = 100, the two modes still converge and give very 
similar results, but they are now at the limit of convergence. The iftc run is 
significantly faster: 56 vs. 205 sec, an almost fourfold improvement. For At — 
200, the non vPtc run does not converge. However, its final iterate temperatures 
still look physical; T r is 0.5% uniformly higher than the corresponding converged 
^tc profile. Our \I>tc implementation has its own limit. The Fig. 5 result, 
where At — 1000, also fails to converge. Nonetheless, the result is physical and 
conserves energy to nearly 11 decimal digits. Non-convergence is evidenced by 
small dips in the matter temperature T m at the cells abutting the left and right 
boundaries. At the end points, T m changes very slowly from one iteration to 
the next. The iterations effectively stall. Although the residuals continue to 
decrease, they have such a slow decay that the run halts when it reaches the 
iteration limit. The run without Wtc and At = 1000 diverges due to negative 
internal energies. To summarize, \E'tc not only decreases the runtime but also 
brings an extra degree of robustness. 

5.3 Expansion of a hot aluminum sphere 

In our opinion, the hardest aspect of code development is integrating a module 
into a multi-physics code and running "real" problems. For us, this implies 
simulations of multiple materials, whose properties are listed in tables, using 
hydrodynamics, heat conduction, radiation modules, and, naturally, AMR. 

For the final test we consider the following problem. An Aluminum (Al) 
sphere of radius 15.5 cm is suspended in air. The initial densities are p = 
2.68118198 and 0.00129 g/cm 3 for Al and air, respectively. Both materials are 
initially at T = 375.936 K. 12 There is initially no radiation energy: £? r |t=0 = 0. 
At t = 0, we inject energy into the radiation field, but only into the domain 
containing Al. The energy is added over 0.1 ns, at which point we have loaded 
a yield Y (erg) into the problem. Energy is added with a Planckian spectrum. 
Unless stated otherwise, the simulations presented in this section use two AMR 
levels; h = 2, 4 cm, and a base grid with h = 8 cm. 



12 Inputs are tailored so that our EOS returns equal pressures for both materials, approxi- 
mately 1 bar. 
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We compare simulations in which radiation transport is modeled by a single 
diffusion equation for the radiation energy density (gray diffusion) to runs where 
the transport is modeled with multigroup diffusion (MGD). We describe results 
where Y w 11 kT and Y w 1 MT. 13 

The problem simulates a strong explosion in air; the parameter choice corre- 
sponds to a nuclear source. The effects are well-known: Zel'dovich and Raizer 
[26] Ch. IX, Brode [5], Landshoff [9]. Initially, radiation dominates the dynam- 
ics: a fast thermal wave propagates through the surrounding air. When the 
wave slows to sonic speeds (of the hot air) , the steep pressure gradient gives rise 
to a strong shock. Finally, hydrodynamics dominates. Salient effects arc similar 
to the simulation of a point explosion using hydrodynamics and nonlinear heat 
conduction (Shestakov [21], "Non-Sclf-Similar-Problem" section). 

Before presenting our results, we summarize them. For the lower yield, gray 
and MGD simulations are very similar. However, for Y = 1 MT, the gray and 
MGD simulations differ significantly and this, we feel, is a new result. Although 
it contradicts established theory (Brode [5]) we believe it to be correct since it is 
explained by examining spectra of the radiation field (see below). Furthermore, 
our MGD result is corroborated by the trusted computer code LASNEX [27]. 

Figures 6, 7 and 8 display densities, temperatures and velocities, respectively. 
Each figure contains three curves. Two are from gray and MGD simulations with 

Y = 11 kT. The third curve is from a simulation using gray diffusion and a yield 

Y = 1 MT. The 1 MT curves are drawn after implementing Sachs scaling, i.e., 
by scaling time and radii by the cube root of the yield ratio Ry = (Yi/Y 2 ) 1 ^ 3 , 
where Y\ — 11 and Y 2 = 1000. Hence, while the Y — 11 results are taken at 
t = 1 ms, the 1 MT results are at t = 4.48 ms and the 1 MT radii have been 
divided by Ry. Figure 6 displays log 1Q (p/p ), where p — 0.00129 g/cm 3 is the 
ambient air density. Although the close agreement displayed in Figs. 6, 7 and 8 
may not surprise, it is indeed remarkable how well the gray scaled 1 MT curves 
compare with the lower yield results. The similarity of the Y = 11 kT gray and 
MGD curves indicates that gray diffusion is adequate for small Y. 

The results in Figs. 6, 7 and 8 are characteristic of an event transitioning from 
a radiation dominated regime to one dominated by hydrodynamics. Figures 6 
and 7 depict a strong shock at r = 31 m separating from a fireball of radius 
26-27 m. 

In order to validate our gray Y = 1 MT simulation, we continue the run 
to t = 7 ms and find good qualitative agreement when we compare with Brode 
[5]. Quantitatively, at t = 7 ms, we find a strong shock at r = 164 m, whereas 
Brode finds it at r as 190 m. Both simulations show a nearly tenfold density rise 
at the shock, while inside the fireball, p 5 • 10~ 5 cm 3 . For the central (r = 0) 
temperature we have T — 2.04 • 10 5 K at t = 7 ms vs. a 2 • 10 5 K for Brode. 
Our fireball radius is 138 m (« 160 for Brode), and our shock temperature is 
1.65 • 10 4 K (a 1.6-1.7 • 10 4 K for Brode). 

We now compare the gray and MGD results for Y = 11 kT yield at the earlier 

13 Using the conversion 4.18 • 10 19 erg/kT, the actual yields arc 10.9731 kT, 0.9870682 MT, 
10.9665 kT, and 0.9862604 MT for the two gray and two MGD runs, respectively. 
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Figure 6: Hot sphere problem. Log of normalized densities; Y = 11 kT yield, 
t = 1 ms, gray and multigroup diffusion; Y = 1 MT gray curve is scaled. 




Figure 7: Hot sphere problem. Matter temperatures T;Y=\\ kT yield, t — 1 
ms, gray and multigroup diffusion; Y = 1 MT gray curve is scaled. 
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Velocities 11 kT, time = 1 ms; 1 MT scaled 
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Figure 8: Hot sphere problem. Velocities; Y — 11 kT yield, t — 1 ms, gray and 
multigroup diffusion; Y = 1 MT gray curve is scaled. 



time, t = 1 fis, when the solution is dominated by radiation. At this time, since 
the thermal wave is supersonic, it suffices to only examine the temperatures T 
and T r , where, for both gray and MGD simulations, T r = (E r ja) x ^ and a is the 
radiation constant. (Although the Al ball has ballooned to nearly 1 m, which 
launches a strong shock at the Al/air interface, there is little separation between 
the interface and the shock. Thus, beyond 1 m, the air density is nearly the 
same as it was initially.) Figure 9, which displays the temperatures, shows little 
difference between gray and MGD. Both models display a fireball extending to 
r = 8.1-8.4 m and a central T ss 2.5 • 10 6 K; both also display the start of the 
shock at the Al/air interface, as evidenced by the spike at r w 0.8 m. 

However, for high yield, the gray and MGD simulations differ dramatically. 
Figure 10 displays T and T r for Y = 1 MT at t = Ifis. 14 We see that for gray 
diffusion, T = T r ; just as for Y = 11 kT. The gray diffusion thermal wave, 
which is still supersonic, has a front at r w 30 m. However, the MGD result 
is strikingly different. Multigroup diffusion lowers the central temperatures by 
more than 10%. More surprisingly, for MGD, T and T r are tightly coupled only 
out to r w 20 m. Beyond that, at T w 8.5 • 10 5 K, T and T r decouple. The 
radiation temperature extends to r w 300 m, which is the free-streaming limit. 

To examine why the high yield gray and MGD simulations differ, we turn 

14 The spatial scale of Figure 10 cannot resolve the small, but nevertheless significant hydro- 
dynamic effects which expand the Al sphere to r ss 80 cm. For MGD, the temperature is not 
monotone w.r.t. to r near the origin. It falls from a central value of 2.6 • 10 6 deg to 1.8 • 10 6 
at the edge of the sphere (due to the rarefying Al) then rises to 2.06 • 10 6 in the air. 



5 SIMULATIONS 



35 



Temperatures T & Tr, 11 kT, time = 1.6-6 S 




Fi gurc 9: Hot sphere problem. Gray and multigroup temperatures T and T^, 
Y = 11 kT, t = lfjjs. 
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Figure 10: Hot sphere problem. Gray and multigroup temperatures T and T r . 
Y = 1 MT, t = l/is. 
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u(x,v) vs. v for various x. t - 1 us 

1 1S , 




Frequency v (keV) 

Figure 11: Hot sphere problem. Spectral radiation energy (erg/cc/keV) vs. 
frequency (keV) at various radii; multigroup physics only; Y = 1 MT, t = 1/is. 



off hydrodynamics and heat conduction, repeat the simulation, and find tem- 
peratures similar to Fig. 10. This is not surprising since the dynamics are 
radiation-dominated. To gain more insight, we examine spectra. Figure 11 dis- 
plays the spectral radiation energy vs. frequency at 5-160 m. Evidently, the 
frequency-dependent air opacity is responsible. High frequency (30-200 keV) 
photons travel largely unimpeded whereas near the origin, the spectrum devel- 
ops a hole at 10 keV. Moving away from the center, the hole progresses to lower 
frequencies so that at 100-200 m, the spectrum consists of two peaks, one at 
the high frequencies, another near the visible range. Since the latter contains 
little energy, the protruding radiation "tongue" of Fig. 10 is due to the high 
frequencies. 

We believe that the difference between the Y = 11 kT and Y = 1 MT MGD 
simulations is due to the factor of 100 between the yields. Because the energy is 
added with a Planckian spectrum, the initial maximum temperatures differ by 
roughly the fourth root, or approximately 3. Since the initial temperatures are of 
order 3-5 keV, the high frequencies have a nearly Wien distribution, v z e~ v l T . 
Hence, we expect the Y = 11 kT spectrum to be e~ v / T /e~ v ' 3T or e - 2 "/ 3T 
times smaller than the high yield case. Substituting T = 3 and v — 100 keV 
gives a very small number. The conclusion is that the Y = 11 kT case has an 
insignificant number of those energetic photons that are not absorbed by air. 

We conclude the section by comparing results of the ID spherical and 3D 
Cartesian versions of our code. We return to running with full functionality, 
i.e., with hydrodynamics, heat conduction, as well as with two AMR levels. For 



6 CONCLUSION/SUMMARY 



37 



1D and 3D temperatures, multigroup, ( = 7.1 ns, V= 11 kT 
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distance (cm) 

Figure 12: Hot sphere problem. Comparison of T and T r for Cartesian and 
spherical multigroup runs; y = 11 kT, t = 7.1 ns. 



the Cartesian simulation, the Al "sphere" is a cube 31 cm per side (in contrast 
to the ID, 31 cm diameter ball.) The difference in volumes implies that the 
initial central, Cartesian temperatures are necessarily smaller in order to have 
the same yield. Figure 12 displays the radial ID results and a x axis lineout of 
the Cartesian run. The agreement of the profiles is self-evident. 

To summarize, we have simulated real-life problems, viz., air bursts with 
yields Y = 11 kT and 1 MT. We've shown that for low Y, gray and MGD 
give similar results. However, for large Y, they differ for early times when the 
dynamics arc dominated by radiation. Our high yield MGD simulation contra- 
dicts results of Brode [5]. However, Brode's pioneering simulations were done 
many years ago when the relatively limited computational resources precluded 
using sophisticated modules such as MGD. 

6 Conclusion / Summary 

We have described a numerical scheme to solve the radiation multigroup diffu- 
sion equations. The scheme is implemented in a radiation- hydrodynamic code 
with the patch-based AMR methodology originally proposed by Berger and 
Oliger [3] for hyperbolic partial differential equations. Our scheme consists of 
two parts. The first, described in Sections 2 and 3, is applied on a level of the 
AMR grid layout and may be adapted to any code. This part consists of adding 
^tc to the "fully-implicit" iterative scheme of Axelrod et al [2]. \I/tc brings 



REFERENCES 



38 



an extra degree of robustness and enhances convergence of the Axclrod scheme. 
We have developed lemmas that determine the minimum magnitude for the iPtc 
parameter r to ensure that the iterations converge and the result is physically 
meaningful. The appropriate magnitude depends on the problem. 

Our implementation of "ftc is not optimal at least for our AMR code archi- 
tecture. In our code, for each AMR level, we compute a single scalar parameter 
r. However, the levels consist of a collection of grids (rectangles in 2D) that 
need not be connected. If the grids are not connected, they form independent 
problems. Hence, it would be more efficient to use different r for disconnected 
grids. 

The second part of our scheme, the sync-solve (SS), addresses a specific need 
of our code, viz., the requirement of having an energy-conserving result on the 
composite grid of multiple AMR levels. For the multigroup equations, this part 
reduces to a coupled system of elliptic equations on the unstructured grid com- 
bining all levels. Since the SS is intended to be a small correction to the result 
of the level solves, we adapted the key element of the "partial temperature" 
scheme of Lund and Wilson [11]. This allowed reducing the multigroup SS to a 
collection of scalar SS's. We were then able to reuse existing software. 

This paper included simulations of three problems. The first two are ideal- 
ized tests of only the multigroup module. The third is a "real" problem, which 
uses the full capability of the code: AMR, multiple materials, etc. The first 
problem was chosen because of its non-triviality and the availability of analytic 
results with which to compare. We obtained excellent agreement and verified 
the convergence properties of the scheme. The second problem illustrated the 
benefits brought by \I/tc. We compared the conventional scheme of Axelrod et 
al [2] with our \Ptc-modified version. For hard problems, \ftc either decreased 
run times or ensured convergence in regimes where the conventional scheme di- 
verged. The third problem showed that our multigroup module has been fully 
integrated into the code and has already extended the scientific frontier. For 
a high yield air burst at STP, we found that gray diffusion gives an incorrect 
result during the radiation-dominated regime because gray fails to capture the 
frequency-dependent effects of the air opacity. 
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Revised S&B table (Bolstad [4]); time t = 1.0, I) = 0.1. Columns 4 and 5 
give maximum, absolute error estimates. Hence, at x = 0, entry T is correct to 
±5.4E-09, i.e., has 8 trustworthy digits. 



B Diagonal dominance; large mean- free- paths 

As noted in the footnote of Section 3.2 (and remarked by a referee), long mean 
free paths may lead to diffusion coefficients that overwhelm the other matrix 
terms. Thus, the estimate for a, obtained in Lemma 2, may be insufficient. The 
matrix diagonal contains three terms of various magnitudes. The first stems 
from the discretization of the d/dt derivative. Because we multiply by At, the 
term equals 1. The second term is due to the coupling coefficient a g ; the term 
equals Atc/lg, where l g = (pKg)^ 1 is the mean free path. The third term is 
the diffusion coefficient, which after including the time step and discretization 
ofd 2 /dx 2 , is of the form Atcl' g /3 h 2 , where l' g is the flux-limitcr-modificd mean 
free path; 

«; = l/[(g- 1 + (3/ l )- 1 (/3 + |A Ug |/n s )], 

where (5 is a small, user-set constant, whose utility will become evident and 
Au g /u g is a normalized difference of adjoining cell-centered values. 

In the limit l g — > oo, the coupling term a g is negligible. Hence, we compare 
the diffusion term with unity. As l g — > oo, l' g no longer depends on l g . After 
factoring a factor of h, the diffusion term is of magnitude, 

(Atc/h)/(f3+\Au g \/u g ). 

The quantity (5~ x plays the role of the maximum number of mean free paths 
allowed, in units of h. If the gradient of u g is not negligible, \Au g \/u g dominates 
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the diffusion term, ff both l g 3> f and \Au g \/u g -C 1, (3 dominates, fn that case, 
we are left comparing 1 to At c/h (5. The parameter j3 is small; we often use 1CP 4 . 
(However, as shown in appendix C, 1CT 4 is too large.) Using c = 3 • 10 10 cm/s 
gives a diffusion term of order 

3-10 14 At/h. (39) 

If this exceeds machine precision, the a estimate of Lemma 2 does not guar- 
antee diagonal dominance. We are now left with problem-specific estimates. 
Clearly, simulations requiring small h or large At are problematic. Luckily, our 
envisioned applications yield reasonable At/h ratios. 

Consider two topics, ICF hohlraums and simulations of the type described 
in Section 5.3. For the former, mesh sizes are rarely less than 0.1 microns, i.e., 
mm(h) = 0(1CP 5 ) cm. Luckily, in ICF, typical total simulations times are of 
order of tens of ns, requiring significantly smaller timesteps. Using max(At) = 
0(1CP 9 ) s, makes (39) of order 10 10 , which, when compared to unity, is six 
orders of magnitude above double precision. 

For applications of the type presented in section 5.3, while timesteps vary 
enormously, so do mesh sizes; hence, the ratio At/h remains moderate. For long- 
time simulations requiring At exceeding 1 s, it is unlikely that it is necessary 
to resolve details less than 100 cm. Substituting these values into (39) leaves 
3 • 10 12 , which, is also resolved by double precision, but just barely. 

C Full physics convergence analysis 

This section presents a spatial convergence analysis of the scheme as it may be 
used in practice. Particular attention is devoted to effects of the flux limiter 
and AMR. In contrast to what was analyzed at the end of section 5.1, here we 
refine about a moving front. 

The exactness of the solution depends upon the magnitudes of At and Ax. 
To ensure that the time step does not dominate the error, we use a very conser- 
vative value for At, which is much smaller than what would be used in practice. 
When At is small, ^tc is not needed. Furthermore, we find that solutions using 
the FI and SI schemes are indistinguishable for our chosen At. We obtain the 
same result by solving nonlinear problems for each time step (FI) as by lineariz- 
ing the equations and solving linear systems (SI). In order to save computer 
time, the simulations in this section use the SI scheme and do not use iff tc. 

To address concerns of a referee, we consider a stringent test and focus 
attention on the problem described in section 5.3. An Al sphere of 15.5 cm radius 
is suspended in air. Initially, both sphere and air are at STP; the radiation field 
is initially zero. We load a Y = 1 MT source (approximately 4 • 10 22 erg) into 
the radiation field only in the region containing Al. The source is loaded into 
a Planckian spectrum over a time interval £ s = 0.1 ns. The interval is so short 
that over its duration the main effects are to raise the radiation field to a high 
temperature and to a lesser extent also increase the matter temperature due 
to coupling. At t = t s , most of the energy is in the radiation field inside the 
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Figure 13: Hot sphere problem. Radiation energy density E r (erg/cc), Y = 1 
MT, t = 0.1 fis, Ax = 0.5 cm. Air region. 



sphere. The radiation temperature T r is largely uniform over the sphere and 
equals approximately 1.3 • 10 8 deg, or nearly 12 keV. 

To highlight effects of the flux limitcr, we examine the solution at t — 10~ 7 s. 
In order to analyze errors due to only our multigroup scheme, we turn off all 
other physics, e.g., hydrodynamics. This yields profiles that are similar to those 
obtained with a "full physics" simulation since at t = 10~ 7 s hydro effects 
should be negligible. (Assuming maximum sound/shock speeds of O(10 7 ) cm/s, 
the most that hydrodynamics can do is push the Al/air interface out a few cm 
while the hot sphere can radiate out to 3000 cm.) 

We examine the total radiation energy density; Fig. 13 displays Er for r > 
40 cm. Inside the Al, E r is much larger than what is shown in Fig. 13; it decays 
steeply from a central value of 2.6 • 10 14 erg/cc, to 1.5 • 10 12 at r = 16, which 
designates the air cell adjoining the sphere. Hence, the radiation temperature 
decays from a central value of T r — 14 • 10 6 to 3.8 • 10 6 deg. The change in slope 
at 1400 cm is explained by examining the temperatures T r and T; see Fig. 14. 
The distance r — 1400 marks the approximate extent of the fireball. However, 
radiation propagates out to r = 3000 cm, then drops sharply; the drop due to 
the flux limiter. 

We identify three distinct regions in the profiles of Figs. 13 and 14. The 
innermost, out to r = 1400 denotes where we can expect the diffusion equa- 
tions to yield an accurate representation of the physics. There, the domain is 
largely optically thick, as evidenced by the close agreement of T and T r . The 
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Figure 14: Hot sphere problem. Radiation T r and matter temperature T (deg), 
Y = 1 MT, t = 0.1 (is, Ax = 0.5 cm. Air region. 



region 1400 < r < 2900 denotes an optically thin region, where the diffusion 
approximation is expected to fail. Lastly, in the region r > 2900 the solution 
depends entirely on a kludge: the flux limiter. However, although at this time 
the limiter is dominant only near r = 3000, it has affected the entire solution 
because the propagation of the front is governed by the limiter and hence all 
cells out to the present position of the front have been traversed by the leading 
edge of the wave. 

For the purposes of the convergence study, we define the results displayed 
in Figs. 13 and 14 as the "converged" solution. We obtain it using a uniform 
grid Ax = 0.5 cm and an initial At = 10~ 16 s. The timestep increases by 5% 
each cycle but is not allowed to exceed 2 • 10~ 12 s. We take 50,183 steps to 
reach the final time. The discretization yields a light-speed Courant number 
C c = cAt/Ax = 0.12. Although the value may seem overly cautious, it is still 
too large. A transport calculation would preclude any signal from propagating 
beyond r = 3016 cm. However, our finest-grid diffusion result yields T r = 
1.8 • 10 6 and 12,000 deg at r = 3016 and 3116 cm, respectively. Although the 
enhanced diffusion of our result may be due to our choice of a limiter, and as 
analyzed by Morel [14] and Olson ct al [15] there are other limiter choices, all 
limiters reduce to discretizing the equation u t = cu x . 

We time-lag the limiter for two reasons. (1) Flux-limiting is a kludge. Thus, 
a time-advanced limiter is not only more complicated to implement but it does 
not yield a more accurate solution. (2) When a front propagates into cold ma- 
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terial, a time-lagged limitcr puts a front slightly behind where a time-advanced 
limiter would place it. 15 So, since the raison d'etre of a limitcr is to retard the 
flow, we time-lag. 

Before describing the convergence study, we take up two topics. The first is a 
truncation error analysis of a diffusion equation with a time-lagged flux limiter. 
If the limiter is dominant, the face-centered diffusion coefficient is D i+1 / 2 = 
[(tt»+i + Ui) /2]/[\ui + i — Ui\/Ax], where i denotes the cell index. After inserting 
this expression into the equation discretized with backward Euler we obtain, 

u' i -u i = (C c /2) [G i+1/2 {u i+1 + - Gi_ 1/2 (ui + ] , (40) 

where Gj+1/2 = {u' i+1 — — Uj| and primes denote the time-advanced 

variable. We ignore the absolute value operator since it only serves to enforce 
flow down the gradient. Since G is of the form f(t+At)/f(t), G»±i/2 = 1+AtF, 
where F has units of inverse time and its leading term is f'(t)/f(t). Expanding 
the LS of (40) about the time-retarded value yields, Atd t u + O(At) 2 , while 
inserting the expansions of G,±i/2 reduces the RS of (40) to 

(CcfflK-Ui+i-Ui-J+AtUt]. 

The term U t is of the form (d t f/f) Ax d x u and has units of u/t. The differ- 
ence (u i+ i — Ui-i) yields 2Axd x u + 0(Ax 3 ). After simplifying, we obtain the 
truncation error of (40), 

d t u + 0{At) = cd x u + 0(Ax 2 ) + C c U t /2 . 

The C c term is important. It shows that even when both At and Ax are 
small, the discretization has an additional error proportional to the light-speed 
Courant number. 

The second topic is related to the parameter introduced at the end of 
section 2. Appendix B shows that for large mean free paths the diffusion co- 
efficient may depend solely on the sum 0/Ax+ |Vu|/u. Consider Fig. 13 and 
the domain 1500 < r < 2500. From values at r = 1500, 2000, 3000, we estimate 
the average of \VE r \/E r to be 3.7 • 10~ 4 . Thus, if a particular is deemed 
sufficiently small for some coarse mesh width, as the mesh is refined, the ratio 
0/ Ax may overwhelm the flux limitcr. The statement has implications for both 
a conventional uniform-grid convergence study as well as one with AMR, since 
we use the same on all levels. Unless stated otherwise, for all runs discussed 
in this section 0/Ax — 10~ 6 . 

We now describe the single-level convergence study. We make five runs 
with successively finer grids, Ax = 8, 4, 2, 1, and for the "converged" result, 
Ax = 0.5 cm. All runs have the same initial At time history. However, max(Ai) 

15 The result may be seen by comparing two face-centered, flux-limited diffusion coefficients 
of the form u/\d x u\ and discretized as (h/2){uo +«i)/|uo — «i |. Assume the front propagates 
into cell 0. Let the time-lagged no = and the time advanced values be u' and u' x with 
u' <C u\. After dividing by h/2, the time-lagged and time-advanced coefficients equal 1 and 
(1 + e)/(l — e), respectively, where e = u /uj. Thus, the time-lagged diffusion is smaller and 
the front does not propagate as far. 
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depends on the size of Ax in order to maintain max(C c ) at 0.12. Hence, max(Ai) 
varies from 3.2 • 1CT 11 to 2 • 1CT 12 between the coarsest and finest runs. 

Our code computes cell-averaged quantities, e.g., E r (erg/cc). On the coars- 
est grid, we obtain 500 values on cells centered at 4, 12, 20, . . ., cm. In order 
to have an equitable comparison, for each run, we also determine 500 average 
values by post-processing. In order to conserve energy, we use volume averaging. 

The error is computed as follows. Let Ei_ \ x denote the averaged energy 
density in the ith coarse cell for a run with mesh size Ax. By defining the 
relative error, 



500 



Ax 



Ei, 0.5) 2 



1/2 



500 



£(£;,o.5) z 



1/2 



(41) 



we obtain 



[e 8 , e 4 , e 2 , ei] = [0.4820, 0.2436, 0.06969, 0.03475] . 

The ratio of successive errors, 

[(e 8 /e 4 ), (e 4 /e 2 ), (e 2 /ei)] = [1.98, 3.50, 2.01]. 

Since the ratios are approximately two, our results suggest first (rather than 
second) order convergence. 

Lastly, we compare results of a run using AMR with an "equivalent" run that 
uses a uniform grid. For the AMR run, we use a base L0 grid with Ax = 8 cm. 
We use two refinement levels, each refines by a factor of two. The run begins with 
At = 10 -16 , and we increase At as before until reaching max(At) = 3.2 • 10~ n . 
Because we refine in both space and time, C c = 0.12 on all levels. The parameter 
= 8 • 10~ 6 ; hence, (3/ Ax varies from 10~ 6 to 4 • 10~ 6 between the coarsest and 
finest levels. The refined levels adapt to the Al/air interface (at r — 15.5 cm) 
and around the position of ma,x[\W(E r )\/E r ]. At the end of the run, the grid 
layout is: < x < 96 and 2944 < x < 3392 for the LI level (Ax = 4) while: 
< x < 48 and 2976 < x < 3360 for L2 (Ax = 2). We compare errors of the 
AMR run as above, by forming cell averages of the 500 cells centered at 4, 12, 
20, . . ., cm. The error on the AMR run is 

e AM R = 0.06963 . 

Since for the AMR run, the finest level Ax = 2 cm, we compare eAMR with 
the error of a uniform-grid run where Ax = 2 cm, i.e., with e 2 = 0.06969. To 
three significant digits, the errors are essentially equal. The uniform-grid run 
uses 2000 cells, while at the end of the simulation, the AMR run has 676 cells. 

To summarize, we find that when the solution depends on the flux limiter, 
spatial convergence reduces to being first, instead of second order. The effect 
reminds us of what happens to hydrodynamic schemes in the presence of shocks: 
in smooth parts of the flow, the solution may be second order convergent, but 
in regions traversed by shocks, the scheme reverts to first order. In closing, we 
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note an additional source of error, which we did not quantify. The simulations 
of this section use real materials whose properties (internal energies, opacities, 
etc.) are given by tables. Table lookups have errors that depend on the schemes 
used to interpolate between table data. 
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